/*
9.27
1.function matrix_calc_taskA,matrix_calc_taskB
2."N" means iteration times
10.2
1.merge the for loop in taskB
3.move memory access variable out of "for loop".
4.add avx256 SIMD 
*/
#ifndef _CALC_H_
#define _CALC_H_

#ifdef __cplusplus
extern "C"{	
#endif

#include<omp.h>//multithread
#include<immintrin.h> //avx
#include<stdio.h>//printf
typedef struct {
	int* 			rowArray;
	const int* 		rowOffset;
	int 			rowArraySize;
	const int* 		columnIndice;
	const double* 	S;
	const double*	valueNormalMatrix;
	double* 		Id;
}TaskMatrixInfoA;

typedef struct {
	const double * 	valueSpiceMatrix;//N*2
	const int* 		rowOffset; 		//n+1
	const int* 		columnIndice; 	//N

	double* 		A;	//N
	double*		 	S;	//nx1
	double*			R;	//nx1
	double* 		H;	//nx1
	const double* 	D;	//N*2
	double* 		IC;	//nx1
	double* 		IG;	//nx1
	double 		alpha;
	int * 			rowArray; //kx1
	int 			rowArraySize;
	void* 			hdl;  
}TaskMatrixInfoB;

void matrix_calc_taskA(TaskMatrixInfoA** listDataList,int N)
{
//int coreNum = omp_get_num_procs();
//omp_set_num_threads(32);
//#pragma omp parallel
//	{
//#pragma omp for
	for (int i = 0 ;i < N; ++i)
	{	
		int row_array_size= listDataList[i]->rowArraySize;	
//		printf("taskA :%d \n",row_array_size);
		TaskMatrixInfoA * list_Data_List=listDataList[i];
		const double *valueNormalMatrix = list_Data_List->valueNormalMatrix;
		const double *S = list_Data_List->S;
		const int *columnIndice = list_Data_List->columnIndice;

		
		for (unsigned int it = 0; it < row_array_size; ++it	)
		{
			const int node = list_Data_List->rowArray[it];
			const int column_num1 = list_Data_List->rowOffset[node];
			const int column_num2 = list_Data_List->rowOffset[node + 1];
			double &id=list_Data_List->Id[node];
			double id2 = id;
			//prepare for avx 
			__m256d m1,m2,m3,sum;			
			sum = _mm256_setzero_pd();
			 
			const int array_left = (column_num2 - column_num1)&0x03;
			int j;
			
			for (j = column_num1;j < column_num2 - array_left;j+=4)
			{
				m1 = _mm256_set_pd(
					valueNormalMatrix[j],
					valueNormalMatrix[j+1],
					valueNormalMatrix[j+2],
					valueNormalMatrix[j+3]
					);
				m2 = _mm256_set_pd(
					S[columnIndice[j]],
					S[columnIndice[j+1]],
					S[columnIndice[j+2]],
					S[columnIndice[j+3]]
					);
				m3 = _mm256_mul_pd(m1,m2);
				sum = _mm256_add_pd(sum,m3);
			}
			for (int j = column_num1;j < column_num2 ;++j){
				id2 +=valueNormalMatrix[j]*S[columnIndice[j]];
			}

			// 多余的数据赋值
			double temp1[4]={0},temp2[4]={0};
			for(int ti=0;ti<array_left;ti++){
				temp1[ti] = valueNormalMatrix[j+ti];
				temp2[ti] = S[columnIndice[j+ti]];	
			}
			//剩余数据乘法
			m1 = _mm256_set_pd(temp1[0],temp1[1],temp1[2],temp1[3]);
			m2 = _mm256_set_pd(temp2[0],temp2[1],temp2[2],temp2[3]);
			m3 = _mm256_mul_pd(m1,m2);
			sum= _mm256_add_pd(sum,m3);

			//取出累加后的数据 
			double tmp[4] __attribute__((aligned(32)));
			_mm256_store_pd(tmp,sum);
			id += tmp[0]+tmp[1]+tmp[2]+tmp[3];
			
//			printf("%.20f\n",id-id2);
			id = id2; 			
		}
	}
//	}
}

//Task B interface
void matrix_calc_taskB(TaskMatrixInfoB** listDataList,int N)
{
//int coreNum = omp_get_num_procs();
//omp_set_num_threads(32);
//#pragma omp parallel
//	{
//#pragma omp for
//	printf("\n\n************************************************\n\n\n");
	for (int i = 0;i<N;i++)
	{
		//指针变量
		TaskMatrixInfoB * list_Data_List = listDataList[i];
		int row_array_size = list_Data_List->rowArraySize;

		const double alpha = list_Data_List->alpha;
		double* A = list_Data_List->A;
		const double* S = list_Data_List->S;
		const int* columnIndice=list_Data_List->columnIndice;
		const double* valueSpiceMatrix=list_Data_List->valueSpiceMatrix;


		// 命名规则：>4循环256d采用大写，<4循环采用小写
		__m256d COND,S_COL,CAP,ALPHA_VEC,TASK_4_5_7_8_IG,TASK_4_5_7_8_IC,TASK_9_A;
		//中间变量
		double 	temp0[4]={0},temp1[4]={0},temp2[4]={0};
		double 	tmp[4] __attribute__((aligned(32)));
		double 	ig;
		double 	ic;
		double 	cond;
		double 	s_col;
		double 	cap;
		int 	p;
		int	u;
		int 	k;
		int	ti;

		ALPHA_VEC = _mm256_set1_pd(alpha);


		for(int it = 0;it < row_array_size;++it)
		{
			//指针变量
			int row = list_Data_List->rowArray[it];
			const int k1 = row << 1;
			double current = list_Data_List->D[k1];
			double charge = list_Data_List->D[k1+1];
			int column_num1 = list_Data_List->rowOffset[row];
			int column_num2 = list_Data_List->rowOffset[row + 1];
									
			double &IG = list_Data_List->IG[row];
			double &IC = list_Data_List->IC[row];
			int array_left= (column_num2 - column_num1)&0x03;


			TASK_4_5_7_8_IG = _mm256_setzero_pd();
			TASK_4_5_7_8_IC = _mm256_setzero_pd();
			TASK_9_A = _mm256_setzero_pd();

			p=column_num1;

			//printf("array_left, total_num: %d %d\n",array_left,column_num2-column_num1);

			for (p = column_num1;p < column_num2-array_left;p+=4)
			{
				k = p << 1;

				// (4)(5) IG矩阵的乘法，即cond矩阵
				COND = _mm256_set_pd(
					valueSpiceMatrix[k+6],
					valueSpiceMatrix[k+4],
					valueSpiceMatrix[k+2],
					valueSpiceMatrix[k]
					);

				// (4)(5) S矩阵
				S_COL = _mm256_set_pd(
					S[columnIndice[p+3]],
					S[columnIndice[p+2]],
					S[columnIndice[p+1]],
					S[columnIndice[p]]
					);
				//IG[row]+=cond*S[col]
				TASK_4_5_7_8_IG = _mm256_add_pd(TASK_4_5_7_8_IG, _mm256_mul_pd(COND,S_COL));

				// (7)(8) IC矩阵的乘法，即cap矩阵
				CAP=_mm256_set_pd(
					valueSpiceMatrix[k+7],
					valueSpiceMatrix[k+5],
					valueSpiceMatrix[k+3],
					valueSpiceMatrix[k+1]
					);		
					
				// IC[row]+=cap*S[col]
				TASK_4_5_7_8_IC = _mm256_add_pd(TASK_4_5_7_8_IC, _mm256_mul_pd(CAP ,S_COL));
				//(9)
				TASK_9_A = _mm256_setzero_pd();
				TASK_9_A = _mm256_add_pd(TASK_9_A,_mm256_mul_pd(ALPHA_VEC,CAP));
				TASK_9_A = _mm256_add_pd(TASK_9_A,COND);

				_mm256_store_pd(tmp,TASK_9_A);
				for (u=0;u<4;u++)
				{
					A[p+u]=tmp[u];
				}
				//printf("1   %f %f %f %f\n",A[p],A[p+1],A[p+2],A[p+3]);
			}

 			if (column_num2-column_num1>=4)
			{
				if (array_left!=0)
				{
					double temp0[4]={0,0,0,0};
					double temp1[4]={0,0,0,0};
					double temp2[4]={0,0,0,0};
					for(ti=0;ti<array_left;ti++){
						temp0[ti] = valueSpiceMatrix[p+p+ti+ti];
						temp1[ti] = valueSpiceMatrix[p+p+ti+ti+1];
						temp2[ti] = S[columnIndice[p+ti]];
					}
					COND  = _mm256_set_pd(temp0[3],temp0[2],temp0[1],temp0[0]);
					CAP   = _mm256_set_pd(temp1[3],temp1[2],temp1[1],temp1[0]);
					S_COL = _mm256_set_pd(temp2[3],temp2[2],temp2[1],temp2[0]);
		
					TASK_4_5_7_8_IG=_mm256_add_pd(TASK_4_5_7_8_IG,_mm256_mul_pd(COND,S_COL));
					TASK_4_5_7_8_IC=_mm256_add_pd(TASK_4_5_7_8_IC,_mm256_mul_pd(CAP ,S_COL));

					TASK_9_A = _mm256_add_pd(TASK_9_A,_mm256_mul_pd(ALPHA_VEC,CAP));
					TASK_9_A = _mm256_add_pd(TASK_9_A,COND);
					_mm256_store_pd(tmp,TASK_9_A);
					for(u=0;u<array_left;u++)
						A[p+u]=tmp[u];

					//printf("2   %f %f %f %f\n",A[p],A[p+1],A[p+2],A[p+3]);
				}

				_mm256_store_pd(tmp,TASK_4_5_7_8_IG);
				current -= tmp[0]+tmp[1]+tmp[2]+tmp[3];

				_mm256_store_pd(tmp,TASK_4_5_7_8_IC);
				charge -= tmp[0]+tmp[1]+tmp[2]+tmp[3];
			}
			else
			{
				for(u=column_num1;u<column_num2;u++)
				{
					k = u << 1;
					ig = valueSpiceMatrix[k]*S[columnIndice[u]];
					ic = valueSpiceMatrix[k+1]*S[columnIndice[u]];
					A[p+u]=valueSpiceMatrix[k]+alpha*valueSpiceMatrix[k+1];
					IG +=ig;
					IC +=ic;
					current -= ig;
					charge -= ic;
				}
				//printf("3   %f %f %f %f\n",A[p],A[p+1],A[p+2],A[p+3]);
			}
			list_Data_List->R[row] = current; 
			list_Data_List->H[row] = charge;
			/********************************************************************/			
			// for(int p = column_num1;p<column_num2;++p)
			// {
			// 	int col = columnIndice[p];
			// 	const int k = p << 1;
			// 	double cond =valueSpiceMatrix[k];
			// 	double cap  =valueSpiceMatrix[k+1];					
			// 	double ig = cond * S[col];
			// 	double ic = cap * S[col];
			// 	//(4)(5)
			// 	IG += ig;
			// 	IC += ic;
			// 	//(7)(8)
			// 	current -= ig;
			// 	charge -= ic;
			// 	//(9)
			// 	A[p] = cond +  cap * alpha;			
			// }
			// list_Data_List->R[row] = current; 
			// list_Data_List->H[row] = charge;
		} 
	}
//	}
}

#ifdef __cplusplus
}
#endif

#endif
